module functions
   use iso_fortran_env, only: dp => real64
   use iso_c_binding
   use petsc_wrapper

contains

   subroutine RHSFunction(u, rhs, t, ctx)
      class(type_domain_vector), intent(in), target :: u
      class(type_domain_vector), intent(inout), target :: rhs
      real(dp), intent(in) :: t
      type(c_ptr), intent(inout) :: ctx

      real(dp), dimension(:), pointer :: up, rhsp
      call u%get_array_read(up)
      call rhs%get_array(rhsp)
      rhsp(1) = -up(2)
      rhsp(2) = up(1)
      call u%restore_array_read(up)
      call rhs%restore_array(rhsp)

   end subroutine RHSFunction

   subroutine ExactSolution(t, u)
      real(dp), intent(in) :: t
      class(type_domain_vector), intent(inout) :: u
      real(dp), pointer :: up(:)
      call u%get_array(up)
      up(1) = cos(t)
      up(2) = sin(t)
      call u%restore_array(up)

   end subroutine ExactSolution

end module functions
program vector_eqn
   use mod_odepack
   use petsc_wrapper
   use iso_fortran_env, only: dp => real64
   use functions

   integer :: i, lrk

   type(type_ode_solver) :: ode_ts
   type(type_domain_vector) :: u, u_exact
   real(dp) :: dt = 0.1_dp, t0 = 0.01_dp
   type(c_ptr) :: ctx
   call petsc_initialize()
   call ode_ts%create("rk_tvd3")
   call u%create_general(2)
   call u%setup()
   call u_exact%duplicate(u)
   call ode_ts%ode_solver%set_solution(u)
   call ode_ts%ode_solver%set_time_step(dt)
   call ode_ts%ode_solver%set_rhs_function(vector_function=RHSFunction)
   call exactSolution(t0, u)
   do i = 1, 99
      do lrk = 1, ode_ts%ode_solver%order
         call ode_ts%ode_solver%substep(lrk, u, ctx)
      enddo
   end do
   call ExactSolution(t0 + (i - 1)*dt, u_exact)
   call u_exact%axpy(-1.0_dp, u)
   call u_exact%print()
   call u%print()
   print *, u_exact%norm()
   call petsc_finalize()
end
